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Abstract. In this paper, we study a parabolic system of three equations which permits to solve 
an optimal replication problem in incomplete markets. We obtain existence and uniqueness of the 
solution in suitable Sobolev spaces and propose a numerical method to compute the optimal strategy. 



I. Introduction 

We study here a parabolic system arising in the resolution of an optimal replication problem in 
incomplete markets. Given a European derivative security with an arbitrary payoff function, the 
optimal replication problem is to find a dynamic portfolio strategy, that is self-financing and comes 
as close as possible to the payoff at maturity date T. In complet markets, such a dynamic-hedging 
strategy exists: the payoff of a European option can be replicated exactly; it is the Black-Scholes 
model (1973) [2]. In pQ , Bertsimas, Kogan and Lo propose a solution approach for this problem in 
incomplete markets.. 

At time r = 0, consider a portfolio of stocks and riskless bonds at a cost Vq and denote by 
0(t), B(t), V(t) the number of shares of the stock held, the value of bonds held and the market 
value of the portfolio at time r. Hence, V{r) — 9(t)P(t) + B(t), < r < T. 
If we note a the volatility and F the payoff function, the value function J is defined by: 

J(r, V, P, a) = min E(((V(T) - F{P(t), a(T))) 2 /(V(r), P(r), a{r))). 

0{s), s>r 

1 /2 

The replication error e(V ) is (J(0, V,P,a)) 1/z and it can be minimized with respect to the initial 
wealth Vq to yield the least-cost optimal-replication strategy and the minimum replication error e* 
is e* =min e(Vo) 

Vo 

In pQ, it has been proved the the value function J is quadratic in V : J = a(V — b) 2 + c and the 
coefficients a, b, c satisfy the following system of partial differential equations: 

,„ „. da k 2 a 2 d 2 a . .da o,o<r 2 f da\ 2 „ 9 . . 



Or 2 da 2 * v 'da r a \da 



db k 2 a 2 d 2 b a 2 P 2 d 2 b , „ „ d 2 b , , db 



dr 2 da 2 2 dP 2 r dadP * v 'da 



Key words and phrases, mathematical model, nonlinear parabolic problem, finite elements method. 
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!1.2) _ (1 _^*<» 

a ao co- 



de k 2 a 2 d 2 c a 2 P 2 d 2 c , 2 _ d 2 c , . dc , 9c 



1.3) -(l-/)itVa 



36 x 2 
&7 



where > 0, p G [—1, +1], ( p is a correlation coefficient) 
gf(cr) = -<5cr(cr -ax) (5 > and g]0, 1[) 

gx(o) = g(a) - 2pkaf(a),g 2 (cr) = g(a) - pkaf(a) 

{—if a < (T 
a fi , p > 0, (p is the drift). 

- if cr > cr 
a 

Remark 1.1. T/ie function f has been modified near m order to be bounded and assure the 
existence of a solution. 

The conditions at the time expiry T are given by: a(T) = 1, b(T) = F(P, a), c(T) = 0. 

Under the optimal replication strategy 6*, the minimum replication error as a function of the 

initial wealth Vq is (J(0))5 = (a(0)(Vo — b(0)) 2 + c(0)) 5 , hence the initial wealth that minimizes 
the replication error is V$ = 6(0) the minimal replication error over all Vq is e* = v/c(0) and the 

least-cost optimal strategy at r = is 9*(0) = ^p(O) + ~j5"q~ (^)* 

Remark 1.2. Exact replication is possible when k 2 (l — p 2 ) = and t/ws corresponds to the following 
cases: 

- Volatility is a deterministic function of time. 

- The Brownian motions driving stocks prices and volatility are perfectly correlated. 



In this paper, we propose a numerical method to compute the solution of equations fll.ip . (jl.2p . (11.31) 
and then obtain the minimal replication error and the least-cost optimal replication strategy. 

To obtain a forward problem, we change the sense of time; we note t = T — r. In order to avoid 
the function a at the denominator, we make the change of unknown ux = ln(a). We also replace a 
by x, P by y, b by u 2 and c by u 3 . 
The preceding system becomes: 
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, . du\ k 2 x 2 d 2 ui . .dui , 2 , 2 1. 2 / dui\ 2 

M) Hi 2~~dx r ~ ^~dx + ~ 2' \~dx' +/W = ° 



<9w 2 ar <9 w 2 x y d u 2 ^ 2 ^ M 2 / \ <9«2 
<9t 2 dx 2 2 (9y 2 cfcrcfy cte 



du 3 k x 6 1 u 3 xydu 3 ^ 2 9 ti 3 , , du 3 . . du 3 
dt 2 <9x 2 2 dy 2 dxdy dx dy 

'dUn 



[1.6) - (1 - p 2 )k 2 x 2 exp(u l ) 



dx 



with the initial conditions: 

Wi(0) = 0; ^2(0) = F(x,y); u 3 (0) = 0.( F is the payoff function). 
The outline of the paper is as follows: 

In section 2, we solve (jl.4p . The different derivative terms will be treated separately in order to 
obtain the L°°- stability of the scheme. We prove the convergence of the numerical solution towards 
a weak solution of the problem. Besides the uniqueness of this weak solution is obtained. 

In sections 3 and 4, we study (|1.5|) . (11.61) . We use a change of unknown which lead to a 
variationnel formulation and obtain the existence of a unique solution in suitable weighted Sobolev 
spaces. These equations are discretized by using a backward Euler method in time and a finite 
element method in space. Numerical results are presented. 

2. Computation of u\ 

2.1. Definition of the numerical solution. In order to solve (jl.4p . we use suitable weighted 
Sobolev spaces, such that no boundary condition is needed in and the function u\ has the correct 
behaviour at infinity. 
To simplify notation, we denote: 

F 1 (x)=f 2 (x), x > 0, A = k 2 {p 2 - i), 

this coefficient may be positive or negative since p is a correlation factor and then p lies in [— 1, +1]. 
The equation ( 11 .411 becomes: 

dui 1 ;2 2 d 2 ui , ,du\ 2 f dui^ ' 
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We will make the following assumptions on the functions / and g\. 
-f G W 1 '°°(R + ), xf G L°°(R + ). 

-F\ is a nonincreasing function, F[ is a bounded variation function; if we denote F the function 
defined by F(x) = xF{(x), x > 0,we get F G L°°(M + ). 

-gi may be written g\{x) = x<p(x) with <p(x) = —5(x — a x ) — 2pkf(x), 5 > 0; so, there exists 
cr 2 > such that is negative and nonincreasing on [cr 2 , +oo[ and bounded on [0, <7 2 ]. 

We define the two constants C\ and c 2 by 

(2.2) ci = sup \g[(x)\, c 2 = sup |0(x)| 

x G [0, cr 2 ] s G [0, o 2] 

We denote by At n the time increment between the levels £„ and t n +i, n > and by vJ[ h the 
approximate solution at the time level t n . This solution will be in a finite-dimensional space Vxh 
which will be defined below. 

The solution u^ 1 at the time level t n+ i is computed in two steps: knowing u^, we compute u\ h 2 , 
approximate solution of 

dui , 9 / <9wi\ 2 

< 2 - 3 > ir +A * (&) =° 

obtained by using an explicit upwind scheme. Then starting with this intermediate value, we use a 
backward Euler method in time to compute u^ 1 ; the second order term in (12.11) is discretized by 
using a Px-finite element method [I] and the linear first order term by an implicit upwind scheme [6] 
in order to get the L°° -stability of the scheme. 

In order to define the finite-dimensional space Vih, we first study the parabolic problem: 

dui 1 , o 9 d 2 ui „ . . 
(2.4) --- kV — + Fl{x)=0 

to obtain a variational formulation in weighted Sobolev spaces. 



2.1.1. Variational formulation of \2.J$ . 
Let us consider the two spaces: 



1 + X 



and 



V xv' 



V 1 = {ve D'(R+)/— — G L 2 (M+), G L 2 (R+] 

1 + x 1 + x 



The space ifi is equipped with the following scalar product: 

r^ 00 y(x*)w(x} 

Vv,w G Hi, (v,w)i= / — zrdx and the associated norm. 

' (1 + x) 







The space Vi with the norm 
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2 /X,\2\ \2 



v 2 (x) x 2 f dv\ 

l"lk=|/ ( TTTxf + (i + ^) 2 W ]<lr 



is a Hilbert space and D(IR + ) is dense in V\ [3]. 
We define on V\ x V% the bilinear form: 

t/ / \ 1 , 2 f + °° dv d ( x 2 \ 
V,, «, G 14, a(t,, «,) = -Jfc | -- ^ J rfx 

or 

, l, 2 f +co x 2 dvdw , , 2 /"+ 00 x cfe , 

a(u,u>) = -« / — ■- — — dx + « / ——wdx. 

2 ,/ (1 + x)^ ax ax J (1 + x) 6 dx 

This bilinear form is continue on V% x V\ and we have the equality 
k 2 f + °° x 2 fdvV , fc 2 /-+ 00 l-2x 



a(t>,t>) = — / - — - — -tt ( ^— ) ax / ^—rV 2 dx, 

K ' 2 Jo il + xf \dx) 2 J (1 + x) 4 

then we get: 



k 2 



Vv G Vi, a(v,v) > '— \\v\\ 2 Vi - k 2 \\v\\* Hl . 

Since the function F\ is in Hi, the following variational problem: 
Find Mi G L 2 (0, T; 14) n C(0, T; Fi) such that: 



(2.5) <( V ^ 

has a unique solution [5]. 



(<1 ,v) + a(ui, v) = —{Fi, v)i, Vv G V\ 



ui(0) = 



2.1.2. Approximation of A2.5\) . 

The finite-dimensional space Vih will be a subspace of Vi defined in the following way: 
Let (xi)o<i<N an increasing sequence (xq = 0). We denote hi = Xi — Xi-i, Jj = (xj-i, Xj), 1 < i < N, 
I N+ i = (x N , +oo) 

Vi h = {v h E C°(R + )/ v h \ h G Pi, 1 < i < N, v h[lN+1 G P } • 



The variable x is the volatility which lies, in practice, in ]0, 1[, so, we may use a constant space step 
h on (0, 1) and an increasing sequence (hi) for x > 1 in order that the number of nodes is not too 
important. 

If v h G Vi h , we denote u» = ^(x;). 

We define on Vih an approximate scalar product: 

K, «, fc ) A = y W + ^ 2 (1 + ^)2^ + ™ ^2(1 + x^ + TT^ J 
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obtained by using the trapezoid method on each interval 7j, 1 < i < N ; the last integral bein£ 
computed exactly. 

We also define the Lagrange interpolate 7Tft.Fi of F\ by: 
n h Fi G Vi h and n h Fi(xi) = Fi(rc f ) = F u , l<i< N. 

The approximate solution of (I2.5p at the time level t n+ i is the solution of : 

+ A^a^™^ 1 , v h ) = (u^ h ,v h ) h - A* n (7r/,Fi,« fc ) ft , Vv/, G Vi ft , 



ti^ = 0. 



This may be written: 



<o +1 = w?o - At„F 10 , 



N 



k 2 T 2 k 2 T 2 



N 



with an = v 1 — , 1 <i < JV-1, ajy 

hi + h i+1 h N + 2(1 + x N ) 

2.1.3. Approximation of the first order terms. We compute now an approximate solution of (12.31) 
by using an explicit upwind scheme. 
Let us denote by the derivative of u™ h 
u n — u n 

We set v£ = 0. 

We shall prove below that the function v% is positive and the function u\ h is negative; we define 
ull 2 G Vi h by: 

(2.6) = < - AAi n x 2 «) 2 •/ A > 0, 



(2.7) ^ 2 = u?, - AAt„x 2 +1 (w" +1 ) if A < 



< z < iV. 

For the linear first order term, since the function gi is not bounded, we use an implicit scheme, 
which will be decentered in order to get a monotone matrix. 
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Finally, the solution G V\h of (12.11) is denned by: 

<o +1 = %o +4 " At„F 10 , 



At 



(2.8) = - At n Fx 



+ ^ («* - In) (u^-u^) 



Ui m 2 - At„Fl 



with 7j = gi(xi), l<i<N, 

c ,0«/ 7! > 0, 
' Mz/7i<0. 

Since g± is negative for x > 02 , we get 5i = 1 if z is large enough. 
The preceding equations may be written by using the derivative v%: 

u\t X + ton ((tti - fc Ti ) < +1 - («i + (1 - 5i)7i) ^m) = %^ " AtnFii, 

< % < N 

2.2. Properties of the scheme. We prove that under a stability condition, the approximate 
solution u\ h is negative and its derivative is positive. 

Let us denote U™ h the vector of M N+1 of components < i < N and F lh the vector of 

components (Fu), < i < N. 

The numerical scheme (12.81) may be written: 

(2.9) (I + to n A h )U^ x = - At n F lh 

where the matrix Ah is tridiagonal and monotone. 



From (I2.6p . (12. 7p and (I2.8p . we get immediately the equations satisfied by v%: 



If v h 2 denotes the derivative of u h 2 , we obtain for 1 < i < N 
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(2.10) 



(2.11) 



and v^ +1 satisfies: 



v n+l + A*n (a _ + a __ i _ ^ + (1 _ 5i _ i)7i _ i} v n+l 



A - (oi-i - 4-i7*-i) <-V - ^ (a* + (1 - W 



hi 



J i+i 



(2.12) 



n+ 



{Fa — Fi,i-x) 



for 1 < i < N 

which may be written: 

(2.13) 



(I + At n B h )V£ +1 = V^~ 2 - At n G lh , 



where is the vector of R of components (uf), 1 < z < iV, -B^ is a tridiagonal matrix (AT x AT), 

Gi/i is the vector of ~R N of components — - — ; — - — , 1 < i < N. 

hi 

Proposition 2.1. If the following stability condition 



(2.14) 



sup 

% > 1 



At n 2 ri 

hi 



< 1 



is satisfied and 

(2.15) At n Ci < 1, 

then the function v% is nonnegative for n > 0. 

Proof: We can rewrite (12.101) as: 



(2.16) 



^ = < ( 1 - A^rfr?) + A^life) 2 « A > 



and we have an analogous formula for A < 0. 
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If (I2.14p is satisfied, we get immediately that > implies v h 2 > 0. Since F x is decreasing, 

the vector V^ +2 — At n Gih is nonnegative and the vector V£ +1 will be nonnegative if / + At n Bh 
is a monotone matrix; this will be true if 1 — ^^(7t — 7«-i) > 0, for 1 < i < N. Since g\ is a 
decreasing function for x > cr 2 , this condition is satisfied for i large enough and if x < o"2„ we get 
1 " ^(7i " 7i-i) > 1 " CiAt n > from d2HHD - 

We deduce immediately the following results: 



Proposition 2.2. Under the hypotheses of proposition \2.1\ the numerical solution u\ h satisfies: 
u n xh < for n > 0. 



Proof: We can rewrite (I2.6P as 



(2.17) ut k = < (l - A^ 2 <) + A^ttfri?^ if \ > 

and we have an analogous equality for A < 0. From proposition ( 12. ip . the function is nonnegative, 
hence if w™ h < 0, we get u^ 2 < 0. 

Since / + At n A h is a monotone matrix and F x > 0, we deduce that u^ 1 < 0. 



Proposition 2.3. Under the hypotheses of proposition [Kl\ the numerical solution satisfies 

II m i^Il°°(1R + ) - tnFl (ty 

for n > 0. 



Proof: We get immediately from (I2.17P : 
matrix, it follows from (I2.9p that 



n+i 
lh ' 



U 



— II w iJIz oo (r+) anc ^ s i nce A-h is a monotone 



\u 



n+l I 



Il°°(m + ) 



< 



lh ' 



U 



L°°(R+) 



+ At n Fi(0) which completes the proof. 



Proposition 2.4. Under the hypotheses of proposition [Ql the function v% satisfies: 

IKII L i (M +) <VMo) 

/or n > . 



Proof: From ( I2.16p . we get if A > 



N 



N 



XAt, 



N 

E 



a? i-l(' u ^-l) ) 
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hence 



n+i 



< vl 



For A < 0, we obtain the same inequality. 



Besides from (J27 



,n+l I 



< 



2D, we get: 



and we deduce the result. 



L!(R+) 



— At n F 1N + At n F 10 . 



We prove now that under some hypothesis on the sequence (hi), the function xv^ is bounded in 
L°°(IR + ) and it is possible to choose the nodes (xi)i<i<Af such that the stability condition is not too 
restrictive. 

We define the function by: 

vi lir v: = Xi v:,i<i<N,vi llN+i = v« N+1 = o. 



Proposition 2.5. If the following stability condition 

(2.18) sup \-^Xi(v™ + «£_!) <lif A > and sup |A| ^x^ + v" +l ) < 1 if A < 



i > 1 



i > 1 



hi 



is satisfied and if the sequence (/i«)i<j<Ar satisfy: There exists a positive constant c such that 

(2.19) 

and 
(2.20) 



1 i 



_Xi hi+x hi ^ _Jm^ x < . < N _ 1 

%i+i hi + h i+ i Xi hi + 

Xn-i h N 



XN 



h N 

- > — c — 

/ijv + 2(1 + x N ) x N h N + hjsr-i ~ x N ' 



then if At < At , At depending on c, c±, C2, the following estimate holds: 



2.21 



\\ v h\\L°°(K , 

for n > and C is a constant depending on c, c\,c^. 



U 



Proof: We get immediately from (I2.10p 

C h = v? - X^x t (Wf - (d) 2 ) if A > 
which may be written: 

«r"'=*r(i-A^(*r+Ci) 

and if (12.181) is satisfied, we obtain: 



+ A^(^+CiR-i 



(2.22) 



< lit) 



L°°(R+) 



h Nl°°(M+) 



If A < 0, we get from fl2~TTD 
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and the estimate (12.221) holds if (12.181) is satisfied. 

Further the following equality results of (I2.12p for 1 < i < N 

+ ^ («i + «i-i - 5,7, + (1 - *i-i)7i-i) ^ n+1 



(2-23) - ^ - 8i-iji-i) -^v^ - ^ (<* + (1 - fc)7i) — 



4 At r 



which may be written 



(j + At n c fe )y fe " +1 = v; n ^ - At n F h 



where F h is the vector of components : Fj = Xj ; — - — , 1 < z < N. Ch is a tridiagonal 

ft* 

matrix (N x N). 

Besides, we have Cu > 0, Cjj < /or z^j, 1 < z < iV. 

We prove that if (12.191) and (12.201) are satisfied, there exists a positive constant c depending on 

N 

C\, c 2 , c such that: 2_, c v — ~^ ^ — * — J? ^- 
We deduce from 02.21) : 

AT 



5^ Cy = j- \ - Qfi-i— 7» + 7i-i + 7i(l - ^i)-^ + 7t-i*i-i— ) > I <i < N -1 

? 1 



; a N - «at_i 7^ + 7^-1 

^iV V ^AT-l X/V-l 

3=1 

Let us denote A- = 



iv 













X. 






xn-i) 


^v V 


-i+7i(l 





1 < i < N- 1, 
A N = — [a N - ajv-r 



and . S 2 = 1 ( - 7i + 7 i_! + Ti (l - — + Ji-iSi-i— ) 2 < / < A. 
In \ Xi+i Xi-\ 
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N 

We get ^Cij = A\ + A 2 . 

3=1 

From (I2~T9D and (12T20D . it follows that A- > -ck 2 , 1 <i < N. 

2 „ x,r , 



Further , we have : A x = — 0(xi)(<5i + 



2 



and for i>2,A { = h (1 - 5») j-<f>(xi) + $i-i<l>{?i-i) , 

We deduce from (12.21) A 2 > — (ci + c 2 ) /or Xj < cr 2 . 

For a; > o"2, the function g\ is negative, so 5i = 1, A 

ft, 

and A\ > 0, since is decreasing. 



A? 



Finally, we obtain c^- > — c, 1 < i < N, with c = cA; 2 + ci + c 2 . 



3=1 



and ii^ +i iil o °(r + ) ^ i _ L* n + Atn II f IL-(m+)) 

The estimate (I2.2ip follows. 

Let us define now a sequence (xj), satisfying (12.191) and (12. 20j) : 

We set: xi = ih, < i < n with n$h = 1, 
then Xi = e 9h Xi_i for i > n + 1 and 9 > 1 

Weget: - > 0, i < n + 1, 

= 0, z > n + 1 



i 



^ r — XN 1 — 7 > — ; then (I2.20p will be satisfied if h N > ■§, that is 

h N + 2(1 + x N ) x N h N + h N _i 2x N 

N = o(\^A) a[XN = o(l\ 



h J \h / 

Xi ( 1 \ At n 

Besides we have — ^ = O ( — ] 1 < ? < iV, and the stability condition may be written — - — < C, 
hi \h ) h 

that is the classical stability condition for hyperbolic problems. 



Proposition 2.6. Under the hypotheses of proposition HOI there exists a positive constant C de- 
pending on T, f, gi such that for t n < T , the following estimate holds: 
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\ V h Wl°°(R+) — 



Proof: For A > 0, we get from (12.101) : 

At 



and by using (12.181) . we obtain: 
For A < 0, we get: 



L°°(R+) 



< \\i] n \\ 

II 6 7iIIl°°(R+) • 



and by using (12. 181) , we obtain: 



n+4 



— \\ v h\\L°°(R+) ( 1 + 2 |A| Ai n ||^|| i0 o( R +) j ■ 



Besides, it follows from ( 12. 13[) that 



,n+l I 



(l- Cl Ai n )||<- ||LOO(R+) 
This concludes the proof. 



< 



+ ^n\\F[\\ t 



Proposition 2.7. Under the hypotheses of proposition \27b\ there a positive constant C depending 
on T, f, gi such that for t n < T , the following estimate holds: 



Var{y 



hi 



< C. 



Proof: For A > 0, we have the equality: 

At 



v i+l 2 — v s = v'^ - r 



Ai r , 



J i+l 



hi+i 



for 1< i < N - 1, and 



— 11 AT 



Air. 



< +1 - ^ + A-^ (a&«) 2 - *a,-iK-i) 2 ) • 



It follows that: 

V i+1 ~ V i 



i+1 -vt){\-^) + ^_M-^ti) 

-2\At n x t (vf +1 + v?)(v? +l - v?) - AAi n /i m « +1 ) 2 - \At n hi{v?) 2 

with = \^x 2 M +1 + <), 1 < i < N - 1, /*» = 0. 

From the stability condition (I2.18p . we get < /x" < 1 and we deduce 

Var < Var(<;R+) (1 + 4AAi n || v^\\ Loo{R+) ) + 2XAt n \K\\ Loo{R+) |K|| L1(M+) . 



14 MARIE-NOELLE LE ROUX 

We have an analogous estimate for A < 0. 



Furthermore, from (12.121) . we obtain the following equalities: 

' < Xi + (1 - oh - Sjru _ 7»+i - li 

hi h i+ i hi + i 



<+i ~ < + + At, 



Wi 1 - 



-At 



ai+i + (1 - 5 i+1 )% +i 



hi 



i+l 



{<2 ~ <i) ~ 



hi 



J i-1 J 



n+7; n+ 



V i+1 



li+1 ~ li _ li - li-1 

hi+i hi 



v? +1 - At r 



Fi,i+i —F\i Fu — Fl,i_i 



hi 



i+l 



hi 



for 1< i < N - 1. 



<+i ~ + ( TW+i ~ < l ) ~ aN -\ lN -H v^ - v- +1 



h y~N+i 



h 



N 



J N U N-1, 



v n+1 2 _ v n+ l _ At lN - lN - 1 4- At FlN ~~ FlN ^ 



h 



N 



TV 



It follows that: 



N-l 

E 

i=l 



i+l 



I - A^ 7m 7 M Ivjft 1 ~v? +1 \ +v n N +1 < Var(v n h +L ,R+) 



N-l 



+At n J2 



li+1 - li li- li-1 



hi+i 



,n+l 



N-l 



+At n 



i=l 



Fu — F\,i-\ Fu — Fi,i_j 



h 



i+l 



hi 



+ At r , 



Fin — F 1 , 



N-l 



h 



N 



It results from (Q: 1 - At n ^ — > 1 - CiAt, 



N-l 



Besides, we have: 



i=i 



K+i 

li+i ~ li li - li-i 



h 



i+l 



hi 



< +1 KVarisk ^ + )H +1 \\ L oo { 



N-l 



and 



-Pi, i+i — Fu Fu — Fn_i 



8=1 



'i+l 



< Var(F[ 



i. ^^) + ^ 
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It follows that: 

(1 - Cl At n ) Var(v n h +1 ;R + ) < Var(v n h +h ;R + ) + CAt n |K +1 || L1(R+) + At n Var(F[, R+) 

and with (12.21) . we get: 

(1 - Cl At n )Var(v]l +1 ; R + ) < Var(v%; E+)(l + 4AAt n ||^|| LOO(R+) ) 

+2\At n H\\ Loo(R+) K|| £1(R+) + CAt n ||< +1 || L1(R+) +At n Var(F{;R+). 
This concludes the proof. 

2.3. Convergence of the scheme and uniqueness of the solution. 

From all these estimates, we can deduce the convergence of the numerical solution to a weak solution 
and we prove that this solution is unique. 

A function u x is called a weak solution of flUD if «i £ C([0,T]); W^ co {R + )), 
x^ £ L7(0,T;L°°(M+)), and 

/ / u\-^-dxdt — / u\(T)(b(T)dx / / *^--^-(x 2 (b)dxdx 

Jo Jr+ dt J R + 2 J J R+ dx dx 

+ f f g\-^-<\>dxdt — A / / x 2 ( — J (f)dxdt = f f Fi4>dxdt 
Jo Jr+ dx J J R+ \dx J J J R + 

for any function <p with a compact support in [0, T] x R + , <fi £ C 1 (0, T x 
Theorem 2.8. Problem \2. 1\) admits at most one weak solution. 

Proof: Let «i and U\ two weak solutions of (12.1 1) . We denote w = U\ — Ui.The function w satisfies: 
dw 1 , 9 9 <9 2 u; . . dw , <9ui\ 

(2 - 24) « +Al &r) = ' 

Let us denote by ip a function in C^R" 1 ") satisfying: 

< ij}(x) < 1, — 1 */ < x < 1, VO 32 ) — if x > 2 and if) decreasing on (1, 2) 

and we define ip u (x) = if)(~), x > 0, v > 0. 

By multiplying (I2.24|) by -g, and integrating on M + , we get: 

(1+x) 

| / w 2 - — — „ dx] + -k 2 [ I tt- I ~ — - — -Tyib u dx + ~ / w 2 „ ^ 1 ^ X \ ( ^-dx 
2dt\J (1 + x) 2 J 2 y V<W (1 + x) 2 ^ 27 (l + a;) 2 cfe 
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2 p+oo 



dw x 2 dil) v , , 2 r 00 <9w 
w- -TT— — dx — k I -z—w- 



dx (1 + x) dx 



x 



o <9x (l + x)' 



Tip u dx 



(2.25) 



2/ "^Iot^J^-Vo ^H^ + ^JaT^ 



ip u dx. 



We estimate now each term of this equality. 



We have: 



+00 



w 



gi{x) dip u 
(1 + x) 2 dx 



dx > if v > cr 2 since g±(x) < and t/£(:e) — 0- 



Besides, we have: ^(a?) = -^(-) and since xff e C(0, T; L°°(M+)), w G C(0, T; L c 
obtain: 



")), we 



2 



+00 



a; 



dipt 



dx (1 + x) dx 



dx 



C 
v 



< — [C depending on ||x^|| LOO(R+) 



and lliwl 



We estimate the second term of the second member of (I2.25P and we get: 



+°° dw x 

w— —%p u dx 



dx (1 + x) 



<k 2 



+00 



w 



(i+xy 



-ip v dx 



1 1 '' +ro fdwY x 2 



\ dx J (1 + x) 



-ip v dx 



, » (dwV x 2 k 2 [ + °° 2 A 



a > 0. 

Further, we have: 



+00 



2 , d ( gi(x) 



W ^Tx\(l+x^ 



dx 



< 



w 



(l + x) 



g[{x){l + x)-2 gi {x) 



l + x 



dx. 



From the hypotheses on the function g±, we get 



g[{x)(l + x) - 2 9l (x) 



l + x 



< 



4>(x) + X(j)'(x) 



and this quantity is bounded. Then , we obtain: 
w ip u — I — I dx < C I w 



dx \(1 + x) 



(x) — x(p'(x)\ 



(l + x) 



rdx. 



It remains to study the last term of (12.251) : 

Since x-^ and x^ e C(0, T; L°°(M + )), we get: 
ax ax 



ax 



ax <9x / ( 1 + x 



x 



;1p u dx 



< c 



dw 



dx 



x 



w 



(1 



X 



;1pvdx 
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< ak 2 I 7—- — szipudx + ^ 7 - f w 2 — . dx. 







dx ) + Aak 2 J (1 + x) 



We deduce from all these estimates: 



1 d ( f + °° 2 i) v ,\ , 2/ l . [ + °° fdw\ 2 x 2 , , 



2dt\J (1 + x) 2 y v 2 7 ,/ v <9xy (1+ x 



If we choose a < |, we obtain by using the Gronwall's lemma: 

y (1+x) 2 " v - - 



+00 w 2 



and if z/ — > + 00, we get / = 

7o (1 + x) 

and we deduce w = and the problem admits at most one solution. 



We prove now the convergence of the numerical solution to this weak solution and thus, we obtain 

the existence of a solution. 

We define the functions UihAt and vx^At by 

lAL n 

V h At =V n h + t ~^K +1 - t)jf), t n < t < t n+1 . 
lAZ n 

Theorem 2.9. Assume that the hypotheses of proposition Iff. 51 are satisfied. The sequence UihAt 
converges uniformly to the weak solution of 112.1}) on any compact of [0, T] x IR + . 

Proof: The functions (uihAt) are uniformly bounded in C(0, T; ty 1,00 (IR + )). 
Further, for R > 0, we get the estimate: 



?/ n+1 — n n 

u lh u lh 



At T . 



L\0,R) 2 



< ^Var(v^;R + )+C( gi ,R) ||< +1 || L i f 



+ \M R \\K\\l™{JL+) \\ V h\\L\R+) + II^iIIli 

with C(gi, R) = sup |gi(x)|. 
x < R 

Thus, the time derivatives of UhAt are uniformly bounded in L°°(0,T; L 1 (0, R)) for any R > 0. 
So, we can extract from the sequence (UhAt) a subsequence, again labeled u h & t which converges 
uniformly on any compact subset of [0, T] x R + to a function u\ [TJ. 
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The functions v^At are uniformly bounded in C(0, T; BV(R + )); 
besides, we have: 



(2.26) 



and 



,n+l n\ 



,n+l n 



H-\0,R)= SU P . 



< v 



n+l n 



H 2 (0,R) 



sup 



<< +1 -<,^> 



H 2 (0,R) 



At r 



H-\0,R) 



< C 



n+l ni n 

h 



At r 



L\0,R) 



< C{R). 



Then, we can extract from (v h/ \ t ) a subsequence, again labeled VhAt which converges to a function 



dui 
dx 



in C(0,T; L\Q)) for any compact Qcl + 0. 



We get easily that u\ is a weak solution. Since this solution is unique, all the sequence is converging 

tO U\. 



So, we have obtained the following result: 

Theorem 2.10. Problem \2. 11) admits a unique weak solution. 



The following figure represents a(0) with a time maturity [3j equal to 1 for different values of p- the 
values of the other parameters are those proposed in [lj: k — 0.4, 6 = 2, o\ = 0.153, \i = 0.7, <7o = 
0.01. 



1.0 




Figure 1. a with T=l 



3. Computation of u 2 



U2 is solution of 
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du 2 k 2 x 2 d 2 u 2 x 2 y 2 d 2 u 2 , 2 °^ u 2 



at 



2 dx 2 



2 dy 2 



— pkx y 



dxdy 



92{X) 



du 2 
dx 



n_ p 2 )k 2 x 2duidu2 
dx dx 



0. 



The initial condition is the payoff function. If the European option is a put, the payoff function is 
given by: F(y) = Max(E — y, 0) if E is the exercise price. 

3.1. Existence and uniqueness of the solution. 

In order to obtain a bilinear form satisfying Garding's inequality, we make a change of unknown. 

We denote: u 2 = e~ ax u 2 , a > 0. 
The function u 2 is solution of: 



du 2 k 2 2 d 2 u 2 x 2 y 2 d 2 u 2 2 d 2 u 2 du 2 ( 2 2 2 2 2 <9«i 

-a; o — — p/cx y-^— ^ ^— I otk x + (felxj + (1 — p J/c x 



2 <y 



at 2~ dx 2 

2 du 2 fa k 2 i 2\ ; 2 2 

x + + (1 — P )a:K x 



dxdy dx 

-pakx z y—^- — I " " x z + ag 2 (x) + (1 — p z )ak z x z ^^- 
dy \ 2 ax 

with the initial condition: ^(O) = e~ ax F(y). 



u 2 = 0, 



<9x 



We define a variational formulation of this problem. 



Let us consider the following space V 2 defined by: 

V 2 = [ v e V'(Q)/ v E L 2 (ft), xv E L 2 (tt), xv x G L 2 (tt), xyv y G L 2 (fi)} 



x 



with = 
This space with the norm: 



v 2 + xV + x 2 f |^ + x 2 y 2 l-l 

is a Hilbert space and X>(f2) is dense inV^ [3] 
Besides, we have the estimates: 
dv 



dv\ \ 2 



\xv\\ T i 



L Z (Q) 



< 2 



xy 



dy 



and the semi-norm: 



ma) m w 



< 2 



dv 
dx 



x- 



L 2 (n) 





dv 


2 


dv 


2 


( 


dx 


+ 

L 2 (C) 


dy 


L2(Q)/ 



is a norm in equivalent to the norm 



lv 2 



L3j 
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We define on V 2 x V 2 the bilinear form b by : 

k 2 f du d , o , , , 1 



\/u,v e V 2 , b(u,v) = y 
,du d 



+pk Jq X dxdy 

—pak / x 2 y—vdxdy — 

Jn oy Jn 



2 

(yv)dxdy — I [ak 2 x 2 + g 2 (x) + (1 — p 2 )k 2 x 2 



du d 

n dxdx^ v)dxdy+ 2j/Tyd-y {y2v)dxdy 

du\ du 
-^-)—vdxdy 
q OX ox 



^-x 2 + ag 2 (x) + (1 - p 2 )ak 2 x 2 ^p- 
2 ox 



uvdxdy. 



We have proved in the preceding section that the function x— — is in C(0, T; L c 

or 

Besides, we have g 2 (x) = —5x(x — a±) — pkxf(x). 

We denote fi(x) = xf(x) and we assume that f\,f[ € 

It is clear that the the bilinear form b is continue on V 2 x V 2 . 



"))■ 



We prove now that it is possible to choose a > such that this bilinear form satisfies Garding's 
inequality. 



Proposition 3.1. // \p\ < 1 and k < |5 ; there exists positive constants C, c, a depending on 
5, <7i, k, p such that 



VveV 2 , b(v,v) > C \\v\\l — c\\v\\ 2 L 2 



(O) 



Proof: We have the equality: 

6 
i=l 

with 

k 2 f dv d 1 f dv d f 

Bl = Yj n d-xd-x {x2v)dxdV ;B2= 2j n ^dy-dy-^^ ** = pk J n X 

B 4 = — f ( ak 2 x 2 + g 2 (x)) + (1 — p 2 )k 2 x 2 ^— \ v^-dxdy; B 5 = —pak [ x 2 y^-vdxdy\ 
Jn\ dx J dx Jn dy 



d~x dy- M ^ 



B« = - 



2 l 2 

a k 



x 2 + ag 2 (x) + (1 - p 2 )ak 2 x 2 °" : ' :i 



in V 2 

It is easily seen, for ^ > 0, 1 < % < 4 

_k 2 

dx T2 



dx 



v 2 dxdy. 



2 



x- 



L 2 (Sl) 



T IHI r2 



L 2 (n) ' 



Bo 



xy 



dv 
dy 



L 2 (Q) 
f 2 dv dv 

B3 = "U (I y T*T» 



\^'L\n)< 



— xv 2 )dxdy and this term is bounded from below: 
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B 3 \ > 



IpI 



dv 
xy Tr- 
oy 



L 2 (n) 



+ k 2 



x 



dv 
dx 



L 2 (n) / 



~\p\[ei IHli2 (n) + ^- |H^2 (n) Y ci > 0, 



B 4 



- f v 2 (2ak 2 x + g'o(x)) dx + (1 — p 2 )k 2 f %-^-x 2 v^-dxdy. 
2 Jn Jn dx dx 



This term is bounded from below: 



B 4 > J (ak 2 -5) xv 2 dxdy-^\p\k\\f[\\ L oc( R ^ IMI^^ 

2 





dv 


• 


dx 



L 2 (Sl) 



k A du± 
x- 

4e, 



dx 



v ^L 2 (n) 



and we obtain 



B A > -e 2 11^11^2^ - (1 - p 2 )e 3 



x 



dv 
dx 



L 2 {VL) 



{ak 2 - 5f 1 



4e 2 



IpI 



- |^| A; ||^||^oo (M+) + ^— 



du i 



x- 



dx 



L°°(Q) y 



^L 2 (n)- 



B 5 >- Y ak\\xv\\ L 2 {Q) 
ak 2 



B e >a[S- 



XV ^L 2 (Q) ~~ J a ( Sa i + ( l ~ P 2 ) k2 



x- 



du\ 



dx 



+) ^L 2 (Q) 



L°°(Q)' 



)xv 2 dxdy 



and we deduce : 



B«> a 6- 



a 



ak' 



e 4 \\xv\ 



(5a, + (1 - p 2 )k 2 



x- 



du\ 



dx 



L°°(VL) 



) 2 + «IpI^II/i|| j L°°(M + ) II^Hl 2 ^)- 



From all these estimates, we get: 
b(v,v)>(l-\p\)^(l-(l + \p\)e 3 ) 



X 



dv 
dx 



2/ +^(i-H) 

L 2 (tt) 2 



dv 

xyj- 
dy 



L 2 (Q) 



+ [aid 



ak 



■) -\-\p\ Y" e ) IHIl 2 (Q) -^MlW.) 



with e = ei + €2 + e 3 + 64 and c is a positive constant depending on the different parameters. 

ak 2 



Let us denote P(a) — a [8 



1 . ak 

\P\ e- 

2 m 2 
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This polynomial admits two positive real roots if k < |5 and e small enough. This last condition is 
generally satisfied in practice (k = 0.4, 5 = 2). So we can choose a > such that P(a) > 0. 
Therefore by choosing e^, 1 < i < 4 small enough, a such that P{o) > and if |p| < 1, we obtain: 
b(v, v) > C\\v\\l — c IMI^2/qv 

We may write problem (13. ip in variational form: 



(3-1) { (^,v) + b(u 2 ,v)=0,\JvEV 2 



Find u 2 E C(0,T; L 2 (£l)) n L 2 (0, T; V 2 ) 

+ b(u 2 ,v) = 0, 
u 2 (0) = e ~ ax F 

Since u 2 (0) G V2, by using proposition 13.11 we deduce the following theorem [5]: 

Theorem 3.2. Under the hypotheses of Proposition [3A\ problem Ii3. 1\) admits a unique solution 



3.2. Numerical solution. 

For the computation of u 2 ,, we shall use a backward Euler method in time and a finite element 
method in space [I]. 

We define a triangulation of Q in the following way: 

- (xi), < i < N is the sequence defined in section (2.1.2). 

- Let (yj), 0<j<M(y = 0) another increasing sequence. 

We denote kj = yj — yj-\ and we assume that the sequence (kj), 1 < j < M is increasing. 
The domain [0,x N ) x [0,y M ] is divided in rectangles ( 

(Vj-uVi), 1 < i < N, l<j<M 
and each rectangle is divided in two triangles by the first diagonal. We denote 7^ the set of these 
triangles. 

We define T 2 , T h 3 , K NM by: 

T 2 = {Kj/Kj = (x N ,+oo[x (%■_!, yj ), j = 1,M} 

T ft 3 = {Ki/Ki = (xi-uXi) x (y M ,+oo),i= 1,N} 

K NM = (%n, +00) x (y M , +00) 

The triangulation is then defined by: T A = U T fe 2 U T fc 3 U K NM . 

We associate to this triangulation the finite-dimensional space V 2 h defined by: 

= (v h e C°(n)/VK e rl v h 1 K e P x , \/K e T 2 , v h j K e P yl \ 

2h ~ 1 T ^ y% 1 k e p ^ y v h 1 k nm e p o J 

Pi is the space of polynomials of degree < 1 in x, y ; P y \ is the space of polynomials of degree < 1 
in y ; P x \ is the space of polynomials of degree < 1 in x ; Pq is the space of constants. 
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Let a > such that P(a) > 0. 

If Vh G V2/1, it is easily seen that Vh = e~ ax Vh G V2. 

The approximate value of u 2 at the time level t n will be in V 2 h- 

If f ft G V2/1, we denote 5/j = e _2m: t> ^ and we define the bilinear form 6 on V 2 h X V2/1 by: Vm^, t>/i G V^, 



-(l-p 2 )k 2 J x 2 -^-^v h dxdy. 



We have the equality: ityj = V-u/j, G and we define the approximate value u^ 1 

of u 2 (t n+ i) as the solution of the following problem: 

< = F h, 

where Fh is the Lagrange interpolate of F in V 2 h, 

(uh, Vh)h is an approximate scalar product in L 2 (fl) and bh an approximation of b, obtained by using 
numerical integration. 

In Fig 2, we present the variation of u 2 (T) in y (or 6(0) in P) for different values of the volatility; 
the parameter p is null; the exercise price E is equal to 1; the others parameters have the same 
values as in Figl. 




Figure 2. b(0) with T=l 
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Fig 3 represents the variation of u 2 (T) in y for different values of p, the volatilty is equal to 0.3. 




Figure 3. b(0) with sigma=0.3 



4. Computation of u 3 

u 3 is solution of 

du 3 k 2 x 2 d 2 u 3 x 2 y 2 d 2 u 3 2 d 2 u 3 , ,du 3 . . du 3 

-W - - —W ~ P y dxdy- " 9{x) ^ " Xf(x)y ^y~ 

2w2 2 , \ ( <9«2 X J 



= (1 — p )k x exp(ui) . 

\ ox 

with the initial condition: ^(0) = 0. 

4.1. Existence and uniqueness of the solution. As for the computation of u 2) we make the 
change of unknown: u 3 = e~ ax u 3 and u 3 is solution of: 

du 3 k 2 x 2 d 2 u 3 x 2 y 2 d 2 u 3 2 d 2 u 3 du 3 . 2 2 , <9u 3 , 2 , 

"ft" ~ I" ~ ~ " pkX V dxdy- ~^ {akX + ~ ~dy~ {pkaX V + Xf(x) ^ 



, a 2 k 2 x 2 \ /a - x 2 

(4.1) — -u 3 I h ag{x) ] = (1 — p 2 )k 2 x 2 exp{u\ — ax) 



dx 



We define the bilinear form c on V 2 x V2 by: 
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w,v g v 2 , c(u,v) = — J n Q^Q^( x2v ) dxd y + 2 j Q ^dy^^ dxdy 

f du d f du 

+pk / x 2 — — (yv)dxdy — / (ak 2 x 2 + g(x))—vdxdy 
J a dxdy J n dx 

2 „ „ 9« , , f ,a 2 k 2 x 2 



n 



(pakx + xf(x))y—vdxdy — / ( h ag(x))uvdxdy. 

dy Jn 2 



It is clear that this bilinear form is continue on V 2 x V 2 and under the same hypotheses as for b, it 
satisfies: Vt> G V2, c(v,v) > C \\v\\ 2 y — c \\v || 2 r 2/ n \ . 

\/ 2 L \\L) 

We denote = (1 - p 2 )k 2 x 2 e Ul , x > . 

Problem (14.11) may be written in variational form: 



(4-2) { l?!£,v)+c{ik,v) = (Ge- a *,v) 



Find u 3 eC 2 (0,T;L 2 {Q)nL 2 {0,T;V 2 ) 
dt 

63(0) = 

If we assume that ^ G C(0, T; L°°(f2)), the second member is continue on V 2 and we get the 
following theorem: 

Theorem 4.1. If \p\ < 1, k < \8 and a such that c satisfies Garding inequality, problem h4-ty 
admits a unique solution. 



4.2. Numerical solution. To compute U3, we use the same method as for the computation of u 2 . 

If Vh G V 2 h, we denote Vh = e~ 2ax Vh and define the bilinear form on V 2 h x V 2 h by: 

k 2 f du h d(x 2 v h ) I f 2 du h d(y 2 v h ) f 2 du h d{yv h ) 

c(u h ,v h ) = — / — dxdy+ - / x — dxdy + pk / x — dxdy 

2 J n dx dx 2 J n dy dy J n dx dy 

~ / g{x)—^-v h dxdy - / xf(x)y—^v h dxdy 
Jn dx J n dx 

and we have the equality: Wuh,Vh G V 2 h, c(u h ,Vh) = c(uh,Vh)- 
The approximate solution u 1 ^ 1 at the time level t n+ i satisfies: 



(uSfc 1 - u% h ,v h ) h + At n c h (u^\vh) = {G,v h ) h , Vv h G V 2h , 
< = 0. 

where Ch is an approximation of c obtained by using numerical integration. 
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Fig4 represents the variation of u 3 (T) in y (or c(0) in P) for different values of a. The parameter 
p is equal to 0. (In the case p — 1 or p — — 1, u 3 is null). The replication error e* is given by: 
e* = v^(0). 



2.0e-004^ 
1 .8e-004- 
1 .6e-004- 
1 .4e-004- 
1 .2e-004- 
' 1.08-004- 
8.0e-005- 
6.0e-005- 
4.0e-005- 
2.06-005- 
O.Oe+000 — 





■ sigma=0.1 




- sigma=0.2 




- sigma=0.5 



0.0 0.2 0.4 0.6 



1.0 1.2 1.4 1.6 1.8 2.0 



Figure 4. c(0) with T=l 



Knowing u 2 , we can compute the least-cost optimal strategy #*(0) = — — (0) H — — (0) = — — (T) + 

oP P oa ay 

pk du 2 , T , 
V dy 

Fig 5 represents the variation of 9* in P for different values of cr; the parameter p is equal to 0. 
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- sigma=0.1 




- sigma=0.2 




- sigma=0.5 



0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 



Figure 5. Optimal strategy 
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